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Abstract. We have performed three-dimensional molecular dynamics simulation of cohesive granular particles under a plane 
shear. From the simulation, we found the existence of three distinct phases in steady states: (I) a uniform shear phase, (II) a 
coexistent phase of a shear band and a gas region and (III) a crystal phase. We also found that the critical line between (II) 
and (III) is approximately represented by £ oc exp(/3yL v ), where £, j6, y, L v are the dissipation rate, an unimportant constant, 
the shear rate and the system size of the velocity gradient direction, respectively. 
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INTRODUCTION 

The interactions among macroscopic granular parti- 
cles, such as sands, are characterized by a repulsive 
and a dissipative ones. The energy dissipation through 
inelastic collisions causes destabilization of a uniform 
state of the system. It is known that there appear clus- 
ters in such systems [1, 2], which may be understood 
by the hydrodynamic equations [3, 4, 5, 6]. When a 
shear is applied to the granular system, high-density 
region, called "shear band", appears [7]. There exist 
many papers to estimate the transport coefficients by 
using kinetic theory [8, 9, 10, 11, 12, 13] and to ana- 
lyze the pattern dynamics by using continuum mechanics 
[14, 15, 16, 17, 18, 19, 20,21,22]. 

On the other hand, fine powders of submicron order, 
such as tonner particles or interstellar dusts, have attrac- 
tive forces such as an electrostatic force [23]. The exis- 
tence of the attractive force causes some new features 
due to the competition of the gas-liquid phase transi- 
tion [24, 25, 26] and the dissipative structure [7]. For in- 
stance, nucleation process near equilibrium is well un- 
derstood [27], but that under a shear is not well un- 
derstood. In this paper, we try to characterize the non- 
equilibrium pattern formation of fine powders under a 
plane shear based on the three dimensional molecular dy- 
namics simulation. 



MODEL AND SETUP 

The system we consider consists of monodisperse 
N(= 10,000) spheres, whose radius is a and mass is m. 
The system size is L x x L y x L z (L = L x = L y ) and we 
choose .x-axis as the shear direction and y-axis as the ve- 



locity gradient direction. We assume that the interaction 
between particles is described by the truncated Lennard- 
Jones potential 

^) = M(Sr-fe) 6 ) CfJ, (1) 

[0 (else) 

with the well depth e, where r ;; - = |r, 7 | is the distance 
between z'-th and j'-th particles and r c is cut-off length (in 
this study we use r c — 3a). For the dissipative force, we 
use 

F™(r u , VlJ ) = -^®{a-\r i} \){v ir r i} )r i} , (2) 

with the dissipation rate £, where V;; is the relative 
velocity vector of z'-th and y'-th particles, 0(r) is the 
step function which is 1 for r > and for otherwise, 
and r is a unit vector proportional to r. We note that £ 
is a dissipation parameter related to the coefficient of 
restitution e, for example, e = 0.998 for £ = 0.1 and 
e = 0.983 for £ = 1 .0 at the temperature T = lAs. Thus, 
we are interested in weakly dissipative situations. This 
weak dissipation is necessary to reach a steady state. 

In general, the boundary effect is important for granu- 
lar systems. It is known that the results strongly depend 
on the boundary condition [28, 29, 30, 31, 32]. In this 
paper, we adopt Lees-Edwards condition [33, 34] to sup- 
press the boundary effect with the aid of SLLOD algo- 
rithm [34, 35]: 

§ - 0) 
at m 

^ = F,-j Pyi e x . (4) 

This model is known as that the system is relaxed to 
the uniform shear state near equilibrium, where r, = 





FIGURE 1. Three typical patterns corresponding to three 
phases in steady states for p = 0.308: (I) uniform shear phase 
(y* = 0.88, = 1.0), (II) coexistent phase of shear band 
and gas (y* = 0.66, £* = 1.0), and (III) crystal phase (y* = 
0.355, C* = 1-0). 



(xi,y i} Zi),Pi = (pxi,Pyi,Pzi) are respectively the position 
and the momentum of 2-th (1 < i <N) particle, V; = r,, 
y is the shear rate and e x is the unit vector in x direction. 
The force acting on i-th particle is given by 



S ( r i/,Vf 



(5) 



In this paper, we choose Tq = 1 .4e as the initial temper- 
ature, which is slightly higher than the critical tempera- 
ture for the equilibrium Lennard-Jones fluid [24, 25]. The 
system is equilibrated in the absence of the shear and the 
dissipation for t = 200(m<7 2 /e) 1 / 2 . Then we apply the 
shear and the dissipation to the system. 



RESULTS 

At first, we choose the average density p = 0.308, 
where p is defined by Na 3 /L 2 L Z . Note that this density 
corresponds to the critical density for Lennard-Jones 
fluid at equilibrium [24, 25, 26]. We fix the system size 
L x = L y = 52(7, L z = 12(7 for this density. We study 
how the steady state depends on parameters such as the 
density, the shear rate and the dissipation rate. We found 
that there exist three distinct steady phases. The typical 
patterns for these phases are drawn in Fig. 1, where 
y* = y(ma 2 /e) 1 / 2 and £* = ^(m<7 2 /e) 1 / 2 . The uniform 
shear phase (I) can be observed for a larger shear rate, 
and the crystal phase (III) is obtained for smaller shear 
rate. We found that the coexistence phase (II) between a 
shear band and the gas region disappears when the shear 




FIGURE 2. The phase boundary lines in the plane of the 
shear rate y and the dissipation rate £ for the density p = 0.308. 




FIGURE 3. The dependence of the average granular temper- 
aure f g on the shear rate y for several dissipation rates £ (for 
p= 0.308). 



rate y and the dissipative rate £ become small. The phase 
boundary lines for p = 0.308 are plotted in Fig. 2. 

Now let us reproduce the boundary lines between two 
phases. For this purpose, we try to characterize the be- 
havior in terms of a simple physical argument based 
on the granular temperature To = (»i/3JV)X)i=i | v,- — f> | 2 
[36, 37], where v = v(y,f) is the velocity field of the y- 
axis. We measure the time averaged granular temperature 
f g , changing the shear rate y. 

The dependence of T„ on the shear rate y for several 
dissipation rates £ is showed in Fig. 3, where f* is de- 
fined by f g * = fg/e. In the case of the phase (I), the 
granular temperature satisfies f g °< y 2 , which is consis- 
tent with the results of dimensional analysis [37, 38]. On 
the contrary, the granular temperature T s becomes zero 
as time goes on, i.e. f g = in the phase (III). 

Here, we investigate the dependence of f g on y in 
the phase (II). A typical configuration of the particles is 
presented in Fig. 4. We can separate the system into two 
regions: (i) the shear band region, and (ii) the gas region. 

To this end, we consider the density profile p(y), 
which is defined as follows: we count the number of 
particles in the region (j — 1)<7 <y < jo, where j is an 
integer, and calculate the local average density pj in that 
region. If p ; - > p, this region is regarded as a part of (i), 
otherwise the region is (ii). We introduce « s b and T s b as 
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FIGURE 4. Typical configuration in pattern (II) (p = 
0.308, y* = 0.66, = 1.0). The system is separated into two 
regions as a function of the local density p (y) and the average 
density p: (i) the shear band region, and (ii) the gas region. 




FIGURE 7. The dependence of the average granular temper- 
ature T g on the shear rate y for several dissipation rate £ (open 
circles) and the results calculated by using Eq. (7) (solid lines). 
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FIGURE 5. The number of particles in the region (ii) n u vs. 
the shear rate y. 




FIGURE 8. Boundary between (II) and (III) for each density 
p, where the solid line indicates £ °= aexp(j3yL v ) with the 
density dependent a. 



the number and the granular temperature in the region 
(i) respectively. n u and T a are the corresponding ones in 
the region (ii). The average granular temperature of the 
system T g is given by 



N — n u n u 

T« = -T sb + —T u . 

s N N 



(6) 



For p = 0.308, the relationship between the shear rate 
y and the number of particles in region (ii) n a is plotted 
in Fig. 5. From Fig. 5, we approximately found n u oc y. 
The relationships between y and T u and between y and 
r s b are plotted in Fig. 6, where we approximately found 
T u ~ const., T s b « y. From Eq. (6), we may write the 
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FIGURE 6. The shear rate y vs. the granular temperature r s b 
in the region (i) and T u in the region (ii). 



relation 



N-viu 



t s i+^t: 



N N 

*2 



= af -bf 1 + c, 



(7) 



where T* = Tje, T* h = T sh /e and T* = T u /e. Here a,b,c 
are constants with respect to the shear rate while they 
depend on the density and the dissipation rate. 

The dependence of T s on the shear rate y calculated 
from the definition and that from Eq. (7) are plotted in 
Fig. 7. We found that the expression (7) well reproduces 
the results of our simulation. 

Now let us study the dependence of the phase bound- 
ary lines among the phases on the density. For this pur- 
pose, we investigate the behavior of p = 0.108, 0.161, 
0.231, 0.362, 0.430, 0.521 (L Y = 88a, 72a, 60a, 48a 
44a, 40a) in the plane of the shear rate y and the dis- 
sipation rate £ for several p. From our simulation, we 
found that the phase boundary between (II) and (III) is 
describes by 

C*-«exp(j3fL;), (8) 

where L* = L y /a. Here, a and j3 are constants with 
respect to the shear rate and the dissipation rate, but a 



solely depends on the density a = a{p) as 

ce(p) = 453.66p 5 -721.97p 4 + 425.66p 3 

-114.02p 2 +13.62p -0.3746. (9) 

Note that j3 is independent of the density. The physi- 
cal reason of Eq. (8) can be understood as follows. The 
boundary between phase (II) and phase (III) might be de- 
termined by the condition whether the particle trapped in 
the potential can pop out. The energy gained by a shear 
8E is proportional to y 2 Ly, and the granular temperature 
T s b in the shear band is proportional to yL y (the inverse 
temperature is proportional to (yLy)^ 1 ). From these re- 
lations, the probability of pop out p satisfies the rela- 
tion p <x e~Pv L y, which means 1/t~ £t*e-Mi, where 
T is characteristic time scale of pop out and (£t*)~ l is 
that of the dissipation. Therefore this relation becomes 
C oc exp(/37Ly). 

CONCLUSION 

We have clarified the relations between the shear rate 
and the dissipation rate in the system for cohesive par- 
ticles under a plane shear. In the phase (II), we found 
that the dependence of the granular temperature on the 
shear rate is approximately given by T g = ay — by 2 + c. 
We also found that the boundary between (II) and (III) is 
approximately described by £ °< exp(j3 yLy). 

We suppose that the mechanism of the boundary be- 
tween (I) and (II) is related to the stability criterion for 
granular gas [14, 15, 16, 17, 18, 19, 20, 21, 22]. This will 
be a future subject of our study. 
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